import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial.distance import pdist, squareform
from sklearn.neighbors import NearestNeighbors

class RKNN_DPC:
    def __init__(self, k=7, m=None, R=None):

        self.k = k # K近邻数量
        self.m = m # 类簇中心数量
        self.R = R # 分配阈值
        self.centers = None
        self.labels = None
        self.dist_matrix = None
        self.local_density = None
        self.delta = None
        self.gamma = None
        self.X_normalized = None
        self.knn_indices = None
        
    # 执行RKNN_DPC聚类
    def rknn_dpc(self, X):
        # X: 输入数据集

        # 算法1: 确定类簇中心
        self.centers = self._determine_cluster_centers(X)
        
        # 算法2: 分配非中心点
        self._assign_non_center_points(X)
        
        return self

    # 算法1: 确定类簇中心
    def _determine_cluster_centers(self, X):
        """
        步骤:
        1. 数据归一化，计算样本的距离矩阵
        2. 根据公式(6)计算每个样本的K近邻密度
        3. 根据公式(7)和公式(8)计算每个样本的代表点和代表值
        4. 根据公式(9)和公式(3)计算每个样本的局部密度和相对距离
        5. 通过决策图或设定阈值确定类簇中心
        """
        n_samples = X.shape[0]
        
        # 1. 数据归一化 - 每个特征缩放到[0,1]区间
        if X.shape[1] > 1:  # 多于1个特征时才进行归一化
            min_vals = X.min(axis=0)
            max_vals = X.max(axis=0)
            # 避免除零错误
            range_vals = np.maximum(max_vals - min_vals, 1e-10)
            self.X_normalized = (X - min_vals) / range_vals
        else:
            self.X_normalized = X.copy()
        
        # 2. 计算样本间欧氏距离矩阵
        self.dist_matrix = squareform(pdist(self.X_normalized, metric='euclidean'))
        
        # 3. 获取每个样本的K近邻
        nbrs = NearestNeighbors(n_neighbors=self.k+1, algorithm='auto').fit(self.X_normalized)
        distances, indices = nbrs.kneighbors(self.X_normalized)
        self.knn_indices = indices[:, 1:self.k+1]  # 排除自身
        knn_distances = distances[:, 1:self.k+1]  # K近邻距离（排除自身）
        
        # 4. 公式(6)计算K近邻密度
        knn_density = np.zeros(n_samples)
        for i in range(n_samples):
            # 计算K近邻距离之和
            knn_dist_sum = np.sum(knn_distances[i])
            # rho_i = exp(-sum_j(d_ij))
            knn_density[i] = np.exp(-knn_dist_sum)
        
        # 5. 计算每个样本的代表点 (公式7)
        rep_points = np.zeros(n_samples, dtype=int)
        for i in range(n_samples):
            # 找到K近邻中密度最大的点作为代表点
            knn_idx = self.knn_indices[i]
            rep_points[i] = knn_idx[np.argmax(knn_density[knn_idx])]
        
        # 6. 计算代表值 (公式8) - 统计有多少个点以当前点为代表点
        rep_value = np.zeros(n_samples)
        for i in range(n_samples):
            rep_value[i] = np.sum(rep_points == i)
        
        # 7. 计算局部密度 (公式9) - 密度加上代表值
        self.local_density = knn_density + rep_value
        
        # 8. 计算相对距离 (公式3) - 点到密度比它高的最近点的距离
        self.delta = np.zeros(n_samples)
        for i in range(n_samples):
            higher_density_points = np.where(self.local_density > self.local_density[i])[0]
            if len(higher_density_points) > 0:
                # 到密度更高的点的最小距离
                self.delta[i] = np.min(self.dist_matrix[i, higher_density_points])
            else:
                # 如果没有密度更高的点，设为最大距离
                self.delta[i] = np.max(self.dist_matrix[i])
        
        # 9. 计算决策值 gamma (公式5) - 密度乘以距离
        self.gamma = self.local_density * self.delta
        
        # 10. 选择类簇中心
        if self.m is None:
            # 使用gamma值自动确定类簇中心
            # 按gamma值降序排序
            sorted_indices = np.argsort(self.gamma)[::-1]
            
            # 计算相邻点gamma值的差值
            gamma_diffs = np.diff(self.gamma[sorted_indices])
            
            # 找到gamma值差异最大的前几个点
            largest_diff_indices = np.argsort(gamma_diffs)[::-1]
            
            # 确定类簇数量：找到gamma值差异较大的位置
            if len(largest_diff_indices) > 0:
                # 找到第一个较大的差异位置
                cutoff_idx = largest_diff_indices[0] + 1
                self.m = cutoff_idx
                print(f"自动确定的类簇数量: {self.m}")
            else:
                # 默认设置为3个类簇
                self.m = 3
                print(f"未能自动确定类簇数量，使用默认值: {self.m}")
                
            # 选择前m个gamma值最大的点作为中心
            return sorted_indices[:self.m]
        else:
            # 如果指定了类簇中心数量，直接选择gamma值最大的m个点
            return np.argsort(self.gamma)[::-1][:self.m]

    # 算法2: 分配非中心点到类簇
    def _assign_non_center_points(self, X):
        """
        步骤:
        1. 将每个中心点c的K近邻点集合的中心点置为c，并给定一个载入值L=K
        2. for x∈Xsort do (Xsort为以局部密度由大到小排序的X)
        3. if num < R then (num为x的K近邻点集已分配的数量)
        4. 将x放入预备队列Xready
        5. else 根据公式(10)计算x的最大类簇权重
        6. L = L - 1
        7. end if
        8. if L = 0 then
        9-17. 处理预备队列逻辑
        18. end if
        19. end for
        20-22. 处理剩余预备队列中的点
        """
        n_samples = X.shape[0]
        m = len(self.centers)
        
        # 初始化标签，-1表示未分配
        self.labels = np.full(n_samples, -1)
        
        # 步骤1：将每个中心点c的K近邻点集合的中心点置为c
        for i, center in enumerate(self.centers):
            self.labels[center] = i  # 标记中心点本身
            # 将中心点的K近邻也标记为该类簇
            for neighbor in self.knn_indices[center]:
                self.labels[neighbor] = i
        
        # 预先计算所有点的K近邻（如果尚未计算）
        if self.knn_indices is None:
            nbrs = NearestNeighbors(n_neighbors=self.k+1, algorithm='auto').fit(self.X_normalized)
            _, indices = nbrs.kneighbors(self.X_normalized)
            self.knn_indices = indices[:, 1:self.k+1]  # 排除自身
        
        # 确保局部密度已计算
        if self.local_density is None:
            print("局部密度尚未计算")
            return
        
        # 步骤2：获取所有点，按照局部密度从大到小排序 (Xsort)
        # 这里应该是所有点，不只是未分配的点
        Xsort = sorted(range(n_samples), key=lambda x: -self.local_density[x])
        
        # 设置参数
        L = self.k  # 初始载入值L=K
        Xready = []  # 预备队列
        
        # 处理Xsort中的每个点
        for x in Xsort:
            # 跳过已分配的点
            if self.labels[x] != -1:
                continue
                
            # 步骤3：计算x的K近邻中已被分配的点数 (num)
            num = sum(1 for j in self.knn_indices[x] if self.labels[j] != -1)
            
            # 步骤3-7：根据num与R的关系进行处理
            if num < self.R:
                # 步骤4：将x放入预备队列Xready
                Xready.append(x)
            else:
                # 步骤5-6：根据公式(10)计算x的最大类簇权重
                weights = self._calculate_cluster_weights(x, m)
                
                if np.max(weights) > 0:  # 如果存在有效权重
                    max_weight_cluster = np.argmax(weights)
                    self.labels[x] = max_weight_cluster  # 分配到权重最大的类簇
                else:
                    # 所有权重为0的情况，放入预备队列
                    Xready.append(x)
                
                # 步骤7：L = L - 1
                L = L - 1
            
            # 步骤8-18：if L = 0 then 处理预备队列
            if L == 0:
                # 步骤10-17：for xr ∈ Xready do
                i = 0
                while i < len(Xready):
                    xr = Xready[i]
                    # 重新计算num
                    num = sum(1 for j in self.knn_indices[xr] if self.labels[j] != -1)
                    
                    # 步骤11-12：if num < R then continue
                    if num < self.R:
                        i += 1
                        continue
                    else:
                        # 步骤14：根据公式(10)计算xr的最大类簇权重
                        weights = self._calculate_cluster_weights(xr, m)
                        
                        if np.max(weights) > 0:
                            max_weight_cluster = np.argmax(weights)
                            self.labels[xr] = max_weight_cluster
                            # 将xr从Xready移除
                            Xready.pop(i)
                            # 步骤15：L = K
                            L = self.k
                            # 不增加i，因为已经移除了一个元素
                        else:
                            i += 1
        
        # 步骤20-22：处理剩余在预备队列中的点
        for xr in Xready:
            # 根据公式(10)计算xr的最大类簇权重
            weights = self._calculate_cluster_weights(xr, m)
            
            if np.max(weights) > 0:
                max_weight_cluster = np.argmax(weights)
                self.labels[xr] = max_weight_cluster
            else:
                # 如果仍然没有有效权重，分配到最近的中心点对应的类簇
                distances_to_centers = [self.dist_matrix[xr, center] for center in self.centers]
                nearest_center_idx = np.argmin(distances_to_centers)
                self.labels[xr] = nearest_center_idx

    # 计算点到各个类簇的权重 (公式10)
    def _calculate_cluster_weights(self, point_idx, num_clusters):
        weights = np.zeros(num_clusters)
        
        for cluster_idx in range(num_clusters):
            # 找出K近邻中属于该类簇的点
            cluster_points_in_knn = [idx for idx in self.knn_indices[point_idx] 
                                   if self.labels[idx] == cluster_idx]
            
            if cluster_points_in_knn:
                # 公式(10)计算权重
                dist_sum = sum(self.dist_matrix[point_idx, idx] for idx in cluster_points_in_knn)
                weights[cluster_idx] = np.exp(-dist_sum / len(cluster_points_in_knn))
                
        return weights

    def predict(self, X):
            """
            预测新数据的类簇标签
            """
            return self.labels